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ABSTRACT 


A  3-D  solid  state  device  simulator  code  "SOSDOR", 
developed  at  DREO,  is  presented.  The  code  uses  a  seven-point  finite 
difference  scheme  to  discretize  Poisson  and  the  continuity 
equations.  The  equations  are  then  solved  using  the  Newton-Raphson 
iteration  method.  Additional  information  pertaining  to  griding, 
carrier  mobility  and  recombination  models  as  well  as  boundary 
conditions  incorporated  into  the  code  are  also  presented.  The 
source  files  of  the  code  and  the  graphical  I/O  interfaces  are  also 
described.  The  code  was  tested  by  simulating  a  PIN  diode  under  0  V 
bias  and  under  a  20V  reverse  bias  condition.  The  simulation  results 
are  in  excellent  agreement  with  the  results  of  simulation  of  the 
same  device  by  the  industry-standard  PADRE  code. 


RESUME 


Un  code  de  simulation  de  composants  semiconducteurs  a 
trois  dimensions  "SOSDOR",  developpe  a  DREO,  est  presente.  Le  code 
utilise  une  methode  de  difference  finie  de  sept  points  pour 
convertir  en  numerique  les  equations  de  Poisson  et  de  continuite. 
Les  equations  sont  alors  resolues  avec  la  methode  d' iteration 
Newton-Raphson.  De  1 ' inf ormation  additionnelle  au  sujet  du 
grillage,  de  la  mobilite  des  porteurs  et  des  modeles  de 
recombinaisons  aussi  bien  que  des  types  de  frontieres  incorpores 
dans  le  code  est  aussi  presentee.  Les  fichiers  de  code  source  et 
1' interface  graphique  entree/ sortie  sont  aussi  deer its.  Le  code  a 
ete  teste  en  simulant  une  diode  PIN  sans  aucune  polarisation,  puis 
soumis  4  une  polarisation  inverse  de  20  V.  Les  resultats  de  la 
simulation  sont  en  excellent  accord  avec  les  resultats  de 
simulation  du  meme  dispositif  par  le  code  de  norme  industrielle 
PADRE . 
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EXECUTIVE  SUMMARY 


In  the  evolving  technology  of  electronic  devices,  device 
modelling  provides  an  investigative  tool  for  understanding  the 
physics  of  semiconductor  devices.  This  is  important  during  the 
design  stage,  where  device  engineers  can  visualize  the  effect  of 
parameter  change  on  the  electrical  characteristics  of  the  device 
much  faster  then  via  prototyping.  Modelling  is  also  a  powerful 
method  of  analyzing  existing  devices  whose  function  is  perturbed  by 
an  external  stimulus  such  as  ionizing  radiation,  heat  or  some  other 
physical  phenomena.  In  this  case,  the  model  helps  in  understanding 
the  changes  in  the  electrical  behaviour  of  the  device  as  a  function 
of  the  change  in  internal  parameters.  Cause  and  effect  between  the 
phenomena  can  be  established  and  the  potential  problems 
circumvented . 

A  3-D  computer  code,  named  "SOSDOR"  (Solid  State  Device 
Simulator) ,  has  been  developed  at  DREO  to  serve  as  an  aid  to  study 
the  effects  of  ionizing  radiation  on  the  solid-state  semiconductor 
devices.  Given  a  specific  geometry  and  doping,  the  code  determines 
parameters  such  as  potentials  or  currents  inside  a  solid-state 
device.  It  can  also  determine  effects  of  applied  voltage  bias, 
change  of  design  or  effect  of  ionizing  radiation  on  these 
parameters.  The  code  has  been  benchmarked  against  a  well  known 
industry-standard  code  called  PADRE  using  a  PIN  diode  as  the  common 
device.  An  excellent  agreement  between  the  outputs  of  the  two 
codes  is  observed.  The  user  can  display  the  results  in  colour  in 
3-D  or  2-D.  At  the  present,  the  code  is  implemented  on  a  IBM  RISC 
6000  Unix  workstation. 
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1 . 0  Introduction 


Computer  simulation  of  solid  state  semiconductor  devices  is 
a  cost-effective  tool  for  evaluating  the  operation  of  these  devices 
in  adverse  conditions  such  as  are  encountered  in  the  space 
environment.  In  these  environments,  the  device  can  be  exposed  to 
anomalous  radiation  fields  or  unusual  temperature  gradients  that 
may  lead  to  malfunction.  Simulation  codes  can  realistically 
predict  the  behaviour  of  devices  under  such  conditions.  They  can 
demonstrate  the  effectiveness  of  design  changes  used  to  circumvent 
anomalous  behaviour.  Simulators  have  become  an  indispensable  tool 
for  quantifying  the  effect  of  process  and  physical  variations  on 
the  electrical  behaviour  of  semiconductor  devices  during  the  design 
process.  Implementation  and  design  of  new  devices  relies  heavily 
on  utilization  of  numerical  simulators  as  they  can  considerably 
reduce  the  time  necessary  to  optimize  the  functional  design  and 
parameters.  It  is,  therefore,  of  no  surprise  that  considerable 
effort  is  being  focused  on  the  development  of  realistic  3-D 
electronic  device  simulators. 

In  this  paper,  a  3-D  computer  code  "SOSDOR"  (Solid  State 
Device  Simulator)  capable  of  simulating  solid-state  electronic 
devices,  is  presented.  The  code  was  developed  at  DREO  for  the 
purpose  of  studying  the  effects  of  ionizing  radiation  on  electronic 
devices.  An  earlier  2-D  solid  state  device  simulator^ ,  developed 
at  DREO,  was  used  to  study  the  effect  of  ionizing  radiation  on  the 
function  of  semiconductor  devices.  Inherent  physical  phenomena 
such  as  electrostatic  fields  and  carriers  transport,  was  simulated 
by  using  an  analog  resistance  network,  also  known  as  the  TLM 
(Transmission  Line  Matrix)  method.  The  code  utilized  an  electronic 
circuit  solver  (SPICE)  to  calculate  the  node  parameters.  The  main 
drawback  of  the  approach  was  that  the  speed  of  simulation  was 
hindered  by  the  speed  by  which  the  SPICE  program  solved  the 
resistance  network.  Spice  also  had  extremely  large  memory 
requirements  due  to  the  method  employed  to  solve  the  node  matrix 
equations.  As  a  consequence,  the  TLM  method  was  unsuitable  for  3-D 
applications,  as  these  require  about  n3/n2  times  as  many  nodes  for 
similar  resolution.  Typically,  for  a  medium  resolution  problem, 
the  value  of  n  (number  of  nodes  in  one  direction)  is  about  50. 
This  motivated  the  development  of  SOSDOR,  a  3-D  numerical  code, 
which  was  implemented  on  an  IBM  RISC  6000  machine. 
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2.  Basic  Equations  in  SOSDOR 
2.1  The  Carrier  Equations 


The  electrical  characteristics  of  many  solid  state  devices 
can  be  specified  by  the  Poisson  (eqn  (1))  and  the  continuity 
equations  (eqn  (2)  and  eqn  (3)). 


V27  =  -q/e  (p  -  n  +  Nb  -  N&) 

(1) 

f  -  +  G- 

(2) 

f  ■  -R°  +  G> 

(3) 

R  and  G  in  the 
generation  rates 
the  electron  and 

continuity  equations  are  the  recombination 
for  the  minority  current  carriers.  Jn  and  J 
hole  current  densities,  respectively,  and 

and 

are 

are 

given  by 

Jn  =  -g|innW  +  kBT\inVn 

(4) 

Jp  =  -g|AppW  -  kbT\ipVp 

(5) 

In  equations  (4)  and  (5)  ,  n  and  nn  are  the  hole  and  electron 
mobilities,  T  is  the  absolute  temperature  and  kB  is  Boltzmann's 
constant.  The  carriers  densities  in  Poisson's  equation  can  be 
expressed  as  functions  of  "quasi-fermi"  (<pn,  <pp)  potentials,  i.e. 

n  =  i^exp  [  ( Y-<pn)  / kBT] 


p  =  j^exp  [  (<pp-Y)  /kBT]  (7) 

Under  the  equilibrium  condition  (np^^)  ,  the  quasi-fermi 
potentials  are  equal  to  zero,  but  become  nonzero  if  an  excess  of 
carriers  is  introduced  or  if  a  bias  voltage  has  been  applied  to  the 
device. 


2.2  The  Mobility  Models. 

In  SOSDOR,  the  electron  and  hole  mobilities  have  doping, 
temperature  and  field  dependent  components  that  have  been  combined 
into  a  relationship  shown  in  equation  (8) . 
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The  doping-dependent  empirical  mobility  model  shown  in  equation  (9) 
below  was  adopted  from  Caughey^2^. 
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V'Nmin 
1  +(N/Nlef)* 


+  l^Mnin 
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The  mobility  parameters  Nref,  jimax,  ^mii>  and  a  were  extracted^2!  from 
experimental  data  obtained  over  a  doping  range  from  1.0el4  cm-3  to 
1.0e21  -3cm.  In  this  model,  the  values  for  the  a  parameter  used 
were  0.76  for  holes  and  0.72  for  electrons.  The  reference  doping 
densities  (Nref  )  were  6.3el6  cm-3  for  holes  and  8.5el6  cm-3  for 
electrons.  The  other  two  parameters,  namely  Mmax  and  ^min'  are 
intercepts  of  a  curve  fitted  into  the  data.  For  holes  Mmax  =  495  cm2 
V-1  s_1  and  Mmin  =  47.7  cm2  V-1  s-1  and,  as  expected,  larger  values 
were  given  for  electrons,  namely  =  1330  cm2  V-i  s_1  and  /zmip  = 

65  cur  V-1  s-1.  The  reduction  of  mobility  with  higher  doping  density 
is  due  to  increased  scattering  of  current  carriers  at  the  impurity 
centres.  Similarly,  increased  lattice  vibration  at  higher 
temperature  have  the  same  qualitative  effect  as  increasing  doping 
density.  For  holes  and  electrons,  fj,T  is  given  by 

Hr  =  2 . 3e9T~2-7 

and 

HT  =  2.1  e9r-2-5 

The  expressions  in  equations  (10)  and  (11)  are  applicable^  in  the 
temperature  range  between  150  to  450  degrees  Kelvin. 

At  the  surface  of  the  device,  the  effect  of  the  electric 
field  on  mobility  was  implemented  into  SOSDOR  using^  the  relation 


(holes) 


( electrons ) 


(10) 


(11) 


V-E  =  V-0(1+XEJ  -°-5 


(12) 


where  E  is  the  electric  field  perpendicular  to  the  surface  of  the 
device  and  x  is  the  field  parameter.  For  electrons  and  holes,  the 
values  of  %  used  were  1.54e-5  cm  V-1  and  5.35e-5  cm  V-1, 
respectively.  Inside  the  semiconductor  bulk  the  field  dependence 
takes  on  the  formf2^ 
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[1  +  J^]P 


V, 


sat 


with  f3  having  integer  values,  1  for  holes  and  2  for  electrons.  The 
carrier  saturation  velocities  vsat  have  been  setf21  to  9.5e6  cms  1 
for  holes  and  to  l.le7  cm  s-1  for  the  electrons.  In  SOSDOR, 
Boltzmann  statistics  is  implied,  and  hence  the  relationship  between 
mobility  and  the  diffusion  coefficient  is  given  by  Einstein's 
relationship  D  =  (kBT/q)/z. 


2.3  Recombination  Models 

Two  recombination  processes  have  been  implemented  into 
SOSDOR,  namely,  the  Shockley-Reed-Hall  and  Auger  processes.  The 
Shockley-Reed-Hall  process  is  caused  by  the  presence  of  trapping 
energy  level  in  the  band  gap  and  the  Auger  process  is  attributed  to 
band-to-band  recombination.  The  total  recombination  rate  is  the 
sum  of  the  two  processes  which  individually  are  given  by 


R 


pn-rii 


SRH 


x  [n+n^xp 


kBT 


]  +  xn  [p+niexp 


kBT 


] 


(14) 


and 

E auger  =  (Cjjn+c^p)  (np  -  Hi )  .  (15) 

where  Ei  =  intrinsic  Fermi  level 
Et  =  trap  energy  levels 
cn  =  Auger  coefficient  for  electrons 
cp  =  Auger  coefficient  for  holes 
n^  =  intrinsic  concentration  with  the  band-gap 
narrowing  effect  included 
Tn  =  doping-dependent  electron  lifetime 
Tp  =  doping-dependent  hole  lifetime 

The  cn  and  c  coefficients  show  a  slight  temperature  dependence151 
in  the  temperature  range  between  77K  and  400K,  which  is,  at  the 
present,  not  included  in  the  code.  For  the  300K  temperature  cn  = 
2 . 8xl0-51cm6  s-1  and  cp  =  9.9xl0“32  cm6  s_1. 

In  addition  to  bulk  recombination,  surface  recombination  has 
also  been  implemented  into  SOSDOR  via161 


■  ft 

>> 


* 
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Rsuxf  ~ 


(16) 


_ pn-nj _ 

(n+n^exp — j—±-)/sn  +  (p+n^exp — *  c  ^SP 

The  parameters  sp  and  sn  are  the  surface  recombination  velocities 
for  holes  and  electrons  respectively.  These  values  depend  on  the 
treatment  of  the  device  surface  during  the  manufacturing  process 
(eg.  etching,  oxidizing) .  Commonly  accepted  values  for  surface 
recombination  velocities  are  approx  l.Oe2  cm  s"1  for  oxidized 
surfaces  and  l.Oe3  cm  s_1  for  other  surface  conditions ( 7 1 . 


2.4  Boundary  Conditions 

At  the  boundaries  of  the  device,  other  than  under  the  gate 
area  or  at  the  contacts,  a  homogeneous  reflective,  also  known  as 
Neumann' st8]  boundary  is  implemented.  This  implies  that  the 
electric  field  at  the  boundary  has  no  normal  component  (equation 
(17))  and  the  net  current  flow  across  the  boundary  are  zero. 

E-n  =  0  (17) 

At  an  ohmic  contact  the  surface  potential  is  fixed  at  a  value 
equivalent  to  the  sum  of  the  applied  potential  V  and  the  space- 
charge  potential  (Dirichlet  boundary) . 


3.  Finite  Difference  Implementation  in  SOSDOR 


SOSDOR  uses  the  finite  difference  method  to  solve  the  basic 
semiconductor  equations  (1-3) .  The  discrete  form  of  the  first  and 
the  second  partial  derivatives  of  a  function  U  are  derived  from  the 
Taylor  series^9! 


U{x+h2 ,  y,  z)  »U(x,y,  z)  +h2  d7(x,y,  z) 


h22  &U{x,y,z)  i  h2 3  Puix.y.z)  ( 
2 !  dx2  3 !  dx3 


(18) 


U(x-hl,y,  z)  ~U(x,y,  z)  -hi  dU(x’ y'  z) 

ax 


hi2  &U{x,y,z)  hi 3  d3U(x,y,z) 


(19) 


2  ! 


dx2 


3  ! 


dx3 


The  solution  for  3U(x,y,z)/3x  can  be  obtained  if  equation  (18)  is 
divided  by  h22  and  equation  (19)  by  hi2  followed  by  subtraction  of 
the  two.  This  will  eliminate  the  second  derivative  terms  in  the 
series  yielding  the  finite  difference  approximation  for  the  first 
derivative  in  x-direction,  i.e. 
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dU(x,  y,  z)  „  hi  U(x+h2,y,  z)  _  h2  U{x-hl,y,  z) 
dx  ~  h2{hl+h2)  hi (hl+h2) 


(20) 


The  error  in  the  above  expression  results  from  the  truncation  of 
the  series  at  the  third  derivative  term,  hence  the  error 
»  (hlh2/3)d3U(x,y,z) /dx3.  The  first  derivatives  in  the  y  and  z 
directions  are  obtained  analogously. 

The  second  derivative  with  respect  to  x,  shown  in  equation 
(21)  is  derived  in  a  similar  fashion,  again  using  elimination  of 
unwanted  differential  terms.  The  results  are  then  added  to  yield 


d2U{x,y,z)  „2U(x+h2,y,z)  ,  2U(x-hl,y,z )  (21) 

dx2  ~  h2(hl+h2)  hl{hl+h2) 

The  total  divergence  operator,  given  in  equation  (22) , 

is  the  standard  seven-point  formula  which  is  graphically  depicted 

v*u  =  y6  JLu ±.u  (22) 

^  Pi  5 

in  Figure  1.  In  Figure  1,  h±  is  simply  a  distance  between  the 
central  node  and  a  branch  node.  In  the  more  general  case,  where 
the  grid  is  nonuniform,  the  expressions  for  p±  explicitly  take  on 
forms 

p1=h1  (hx+h2)  ,  p2=h2  (hx+h2)  ,  p 3=h3  (h3+h4) ,  .  .  .  p6=h6  (h5+h6)  (23) 

The  scaling  factor  t  is  equal  to  6/h2  if  the  grid  is  uniform  and, 
for  the  case  of  a  non-uniform; 

£  =  —2—  +  —2—  +  — fL_  (24) 

hxh2  h3h4  h5h6 


4.  The  Solution  Method 

The  solution  to  Poisson  equation  requires  finding  the 
solution  to  the  matrix  equation  AT  =  g  where  A  is  the  divergence 
operator  matrix  of  size  (Mx-1) (M^-l) (Mz-1)  x  (Mx-1) (My-1) (Mz-1) ,Y 
is  the  1  x  (Mx-1)  (M^-l)  (Mz-1)  column  vector  representing  the  node 
potentials  and  g,  also  a  column  vector,  is  the  R.H.S.  of  equation 
(1)  .  For  3-D  systems,  due  to  the  large  matrix  sizes  involved, 
instead  of  solving  the  matrix  equation  by  a  direct  method  (eg. 
Gaussian  elimination) ,  an  iterative  method  of  solution  is  used. 
SOSDOR  solves  the  device  equations  iteratively.  The  iterative 
technique  has  an  advantage  over  the  direct  method  in  large  systems, 
because  of  its  simplicity  and  ease  of  implementation.  The  main 
disadvantage  of  this  technique,  however,  is  the  possibility  of  a 
relatively  slower  rate  of  convergence.  The  device  equations  are 


solved  by  the  Newton-Raphson  iteration  scheme  for  which  the 
recursion  formula  is  given  by 


Vs+l  = 


F(tys,n,p) 
F'(tys,  n,p) 


(25) 


*  Here,  s  is  the  iteration  loop  number,  F(¥,n,p)  is  the  equation  (1), 
(2)  or  (3)  in  the  implicit  form  and  f' (7 ,n,p)  is  dF(¥ ,n,p) /d¥.  The 
criteria  for  the  sequence  to  converge  has  been  discussed  by 
Vilenkin^101  and  is  not  discussed  here  in  detail.  Some  salient 
points,  however,  should  be  mentioned.  First,  difficulties  can 
arise  if  F*  approaches  zero.  In  the  case  of  the  present 
application,  application  of  equation  (25)  does  not  pose  any  danger 
since  F'  is  an  exponential  function.  Second,  the  number  of 
iterations  required  for  convergence  vary  from  node  to  node  and  from 
case  to  case.  A  very  rapid  convergence  with  as  few  as  five 
iterations  has  been  observed  at  locations  within  the  device,  remote 
from  active  areas  (junctions)  and  as  many  as  300  in  the  active 
areas  of  the  device.  Provided  that  the  first  approximation  is 
close  to  the  solution,  the  method  converges  quadratically  according 
to  j  Ys+1-Y ] < j ¥S-Y ] 2 .  This  became  evident  in  the  neutral  regions  of 
the  device,  where  the  initial  potential  and  carrier  density 
estimates  (obtained  from  doping  densities)  were  close  to  final 
values . 

The  time-dependent  continuity  equations  were  first  converted 
to  Crank-Nicolson  implicit  form  and  then  solved  by  Newton-Raphson 

*  iteration.  Crank-Nicolson  implicit  discretization  is  more  stable 
than  a  simple  explicit  discretization  methods!91.  Attempts  to 
solve  the  continuity  equations  via  explicit  Richardson  method  have 

\  resulted  in  instability  and  failure  after  several  iterations.  The 

undesirable  relative  complexity  of  the  implicit  form,  lies  in  the 
fact  that  the  present  unknown  variable  (time=t+l)  is  expressed  by 
other  present  unknown  variables,  as  well  as  variables  from  the 
previous  time  step  (time=t) .  The  continuity  equation  for  holes, 
for  example  in  the  implicit  discrete  form  is 

(i-i)  s6i=i— Pi+1-5pt+1+is6i=i— p/-^p  fc)  =  (l-D  gt+1  (26) 

Pi  Pi 


where  1  is  a  weighing  factor,  p^p  is  the  concentration  of  holes  at 
the  particular  node,  t  is  the  time  and  g  is  the  R.H.S.  of  equation 
(2).  For  example,  1  is  equivalent  to  0.5  in  case  of  Crank-Nicolson 
method,  and  is  equivalent  zero  for  the  explicit  form. 


5.  The  Computer  Code 

Sosdor  is  a  3-D  computer  code  written  in  C  language  for  the 
IBM  RS-6000  Unix-based  workstation  and  runs  in  an  X-windows 
environment.  The  code  is  composed  of  three  main  modules,  namely 
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Figure  1  Seven-point  3-D  finite  difference  computational  scheme 

i)  Input  module 

ii)  Processing  module 

iii)  Output  module 

In  the  Input  module,  the  user  enters  information  ralated  to 
device  layout,  doping  profiles,  electrical  connections  and  the 
device  computational  grid.  The  module  consists  of  a  graphical 
drawing  package  and  a  set  of  standardized  device  parameter  input 
files  that  must  be  generated  by  the  user  prior  running  the  code. 

The  menu-driven,  graphical  drawing  package  allows  the  user  to  / 

graphically  design  the  computing  grid  inside  a  graphical  window. 

For  a  3-D  grid,  the  user  is  required  to  draw  a  minimum  of  two  2-D 
grids  in  x-y  and  x-z  or  y-z  planes  to  overlay  the  device.  Node 
distances  are  then  automatically  placed  inside  the  grid  description 
file  and  subsequently  read  by  the  Processing  module.  As 
alternative  to  a  graphical  mesh  design,  the  user  also  has  an  option 
to  manually  enter  the  node  distances  in  the  x,  y  and  the  z 
direction  into  the  grid  description  file  called  meshdata.  This 
method  of  grid  generation  was  found  to  be  satisfactory  in  a  lot  of 
applications  (dynamic  re-griding  at  this  point  is  not  available) . 

The  file  named  design  contains  the  information  about  the  number  of 
diffusion  wells  and  the  extent  of  these  wells  in  terms  of  node 
numbers.  The  file  electrical  holds  information  about  the  voltage 
bias  of  the  device.  Detailed  information  about  the  files  and 
subroutines  in  the  graphical  part  of  the  Input  module  is  given  in 
greater  detail  elsewhere t 11 1 . 

S 

The  Processing  module  processes  the  device  input  parameters, 
solves  the  transport  equations  and  puts  the  results  into  a  set  of 
files  called  potentials ,  holes  and  electrons.  The  Processing 
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module  is  composed  of  sub-modules  start ,  device ,  functions ,  grid 
and  variable-related  header  files.  At  the  beginning,  the  user  is 
prompted  to  reply  whether  the  device  to  be  simulated  is  new.  The 
replay  directs  the  code  into  either  the  start  submodule  or  device 
submodule.  If  device  simulated  is  not  new,  the  input  module  need 
not  be  used.  The  Processing  module  will  simply  read  the  existing 
output  data  files  obtained  earlier,  and  will  examine  the  effect  of 
changing  voltage  bias  or  the  effect  of  ionizing  radiation  on  the 
potential  and  the  distribution  of  current  carriers.  This  feature 
is  useful  for  the  cases  where  unbiased  devices  are  stored  and  then 
retrieved  for  further  simulation  study.  The  updated  output  files 
receive  the  prefix  "device".  The  functions  submodule  contains  all 
the  functions  definitions  that  are  called  by  other  submodules  of 
the  Processing  module. 

The  Output  module  reads  the  potential  or  current  carriers 
distribution  file  and  the  grid  description  file  to  generate  a 
coloured  graphical  output.  The  densities  or  magnitudes  are 
indicated  with  21  colour  chart  available  with  displays.  There  are 
several  modes  of  display  are  available, 

i)  3-D  whole  device 

ii)  3-D  surface 

iii)  2-D  map 

iv)  2-D  exploded  view. 

After  entering  the  data  file  names,  the  main  (default) 
display  shows  the  device  of  interest  in  3-D.  The  3-D  surface 
displays  one  layer  of  the  device  at  the  time.  The  3rd  dimension 
relates  to  the  magnitude  of  the  variable.  The  2-D  map  displays  one 
layer  of  the  device  at  the  time.  The  2-D  exploded  view  shows  the 
whole  device  by  showing  all  planes,  separated  by  small  space, 
simultaneously.  To  activate  the  other  3  displays,  the  user  makes 
a  selection  from  the  menu  bar  located  at  the  top  of  the  main 
display. 


6.  Computational  Results  and  Discussion 

The  device  geometry  considered  for  the  simulation  is  shown  in 
Fig.  2a.  The  device  is  a  MRD500  PIN  type  diode^12],  a  planar 
silicon  structure  with  both,  the  p++  and  n  regions  heavily  doped  in 
relation  to  the  I  (intrinsic)  region  of  the  device.  All  the  diode 
regions  were  considered  to  be  uniformly  doped,  the  junctions  were 
assumed  to  be  abrupt  and  Shockley-Read_Hall  type  recombination 
process  was  assumed  throughout  the  structure.  Device  parameters 
used  in  the  simulation  are  given  in  Table  1. 

The  grid  used  in  this  simulation  is  three  dimensional, 
nonuniform  and  orthogonal.  The  grid  has  45  nodes  uniformly 
separated  by  11.5  microns  in  the  x  and  y  directions  and  has  35 
nodes  in  the  z  direction  (vertical)  with  node  distances  varying 


n 


Table  I  Physical  parameters  characterizing  the  PIN  diode 
structure 


p++  I 


Doping  cm-3 

Na=7.88el8 

Nd  =4.18el3 

Nd=4.18el5 

Thickness 

4  /zm 

23  /zm 

25  /zm 

Bias 

p++contact=0V 

n  contact=+20V 

Contacts  A1  506  /zm  x  506  /zm  Au  506  /z m  x  506  /zm 

Die  506  /zm  x  506  /zm  x  52  /zm 

Grid  3-D  nonuniform,  orthogonal  45  x  45  x  35 

Temperature  300  K 


from  0.075  microns  to  10  microns,  the  finer  mesh  being  utilized  in 
the  junctions  areas  (Figure  2b).  A  constant  time  step  of  0.01  ps 
was  used  in  all  calculations.  The  temperature  was  assumed  to  be 
300°  K. 

Figures  4a,  4b  and  4c  show  the  potential  and  carrier 
distributions  under  a  zero-bias  condition  as  a  function  of  device 
depth.  The  built-in  potential  varies  from  -0.4855  V  at  the  p++ 
side  of  the  diode,  to  0.334  V  at  the  n  side  of  the  device.  One 
observes  that  the  potential  change  occurs  almost  exclusively  in  the 
intrinsic  region  (I)  of  the  device.  This  is  to  be  expected,  due  to 
the  low  doping  (high  resistivity)  of  the  I  region  in  relation  to 
P++  and  n  (Table  1) .  The  corresponding  hole  and  electron  density 
profiles  also  change  in  the  I  region  of  the  device.  Holes  have  two 
equilibrium  minority  carrier  concentrations,  one  in  the  I  region 
and  one  in  the  n  region.  Both  satisfy  pn  =  ni2,  the  ^  being  the 
intrinsic  density  of  carriers.  The  doping  density  of  the  n  region 
is  about  two  orders  of  magnitude  higher  the  doping  level  of  the  I, 
however,  the  value  of  n.^  is  essentially  the  same  across  these 
regions  of  the  device.  As  the  doping  level  approaches  approx  le16 
cm-3  there  is  a  prominent  decrease  in  the  energy  band-gap.  This 
results  in  an  increase  in  the  value  of  n.^  concentration.  This 
band-gap  narrowing  effect  is  important  in  the  p++  area  of  the 
device  (Na=7.88el8  cm-3),  where  the  value  of  nL  is  about  5  times 
higher  then  in  both  the  I  and  n  regions.  Figure  3  shows  how  band- 
gap  narrowing  affects  the  value  of  intrinsic  concentrations  of 
carriers  in  regions  of  the  device. 

The  results  of  the  PIN  diode  simulation  by  SOSDOR  were 
compared  with  the  simulation  results  of  the  same  device,  carried 
out  at  Renssealer  Poly.  Inst.  (Troy,  N.Y.)  using  the  well  known 
solid  state  device  simulation  code  PADRE. 


CARRIER  CONCENTRATION  (cm- 3) 


DEPTH  (MICRONS) 

Figure  3  Doping-related  band  gap  narrowing  effect  in  the 
intrinsic  concentration  of  carriers 


Figure  4a  Built-in  potential  of  the  PIN  diode 
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Figure  4c  Electron  distribution  inside  the  PIN  diode 
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The  purpose  of  this  comparison  was  to  validate  the 
correctness  of  the  mathematical  approach  selected  and  its 
implementation  into  the  algorithm  of  the  code.  From  Figure  5,  one 
easily  can  conclude  that  there  is  an  excellent  agreement  between 
the  results  from  the  two  codes. 

Figures  6a,  6b  and  6c  are  3-D  plots  showing  simulation 
results  when  the  diode  is  reverse-biased  at  20  V.  The  regional 
contacts  are  connected  to  the  diode  via  metallization  (Figure  2a) 
located  at  the  top  and  the  bottom  surface  of  the  device.  The 
contacts  have  the  same  area  as  the  die  dimension,  namely  506  x 
506  Mm.  The  contacts  were  considered  to  be  ohmic,  so  the  surface 
nodes  at  these  areas  were  set-up  as  Dirichlet  boundaries.  For  the 
highly-doped  regions  of  the  device  such  as  the  p++  region, 
metallization  is  applied  directly  on  top  of  the  semiconductor 
material.  For  the  n  region,  however,  the  ohmic  contact  would  be 
accomplished,  for  example,  via  metal-n+  method  (not  shown  in  Figure 
2a)  .  Oscillations  in  the  intermediate  results  occurred  when  the 
bias  voltage  was  applied  suddenly,  which  in  some  cases  caused  the 
simulation  results  to  diverge.  The  reason  for  such  an 
instability  was  attributed  to  the  sudden  dominance  of  the  drift 
current  over  the  diffusion  current  inside  the  I  region  of  the 
device.  The  problem  was  remedied  by  ramping  the  bias  voltage  in 
small  increments.  Specifically,  in  this  situation,  0.2  V 
increments  were  found  satisfactory. 

Figure  6a  refers  to  the  potential  profile  while  Figures  6b 
and  6c  relate  to  hole  and  electron  densities  respectively.  Due  to 
the  x-y  symmetry  of  the  device,  the  potential  profile  and  the 
carrier  densities  are  constant  across  the  horizontal  planes,  so 
only  vertical  profiles  are  given.  All  the  potential  change 
occurring  within  the  device  is,  again,  confined  within  the  I  region 
of  diode,  thus  resulting  in  an  area  of  strong  electric  field.  The 
I  region  also  becomes  virtually  depleted  of  mobile  carriers.  In 
this  mode  of  operation,  the  PIN  diode  can  be  utilized  as  an 
ionizing  radiation  detector.  Note,  that  although  the  p++  region  is 
only  4  microns  thick,  the  proximity  of  the  contact  potential  and 
the  very  high  doping  profile  of  this  region  prevents  any 
significant  modulation  of  device  parameters  in  this  region  of  the 
device. 


7.  Conclusion 

The  SOSDOR  code,  developed  at  DREO  for  simulation  of 
semiconductor  devices,  has  been  described.  The  simulation  of  an 
MRD-500  PIN  diode  has  been  compared  with  the  simulation  results 
obtained  with  an  industry  standard  code  PADRE.  The  code  can  be 
utilized  for  design  of  radiation  detectors,  solid  state  device 
characteristics  studies  and  studies  of  radiation  effects  on  device 
function.  In  the  future,  the  code  will  be  extended  to  include  a 
GUI  (graphical  user  interface)  and  a  re-griding  algorithm. 
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DENSITY 


(b) 


DEPTH  (microns) 


(c) 


0  12  24  36  48  60 

DEPTH  (microns) 


DEPTH  (microns) 


Figure  5  Comparison  of  the  simulation  results  between  SOSDOR  and 
PADRE. 

a)  Built-in  potentials,  b)  Hole  density,  c)  Electron 
density 


Figure  6c  Distribution  of  electrons  inside  the  PIN  diode  for 
20V  reverse  bias. 
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